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Abstract 



This paper examines the onset of the viscous overstabihty in dense particulate 
rings. First, we formulate a dense gas kinetic theory that is applicable to the Sat- 
urnian system. Our model is essentially that of Araki and Tremaine (1986), which 
we show can be both simplified and generalised. Second, we put this model to 
work computing the equilibrium properties of dense planetary rings, which we sub- 
sequently compare with the results of A^-body simulations, namely those of Salo 
(1991). Finally, we present the linear stability analyses of these equilibrium states, 
and derive criteria for the onset of viscous overstability in the self-gravitating and 
non-self-gravitating cases. These are framed in terms of particle size, orbital fre- 
quency, optical depth, and the parameters of the collision law. Our results compare 
favourably with the simulations of Salo et al. (2001). The accuracy and practicality 
of the continuum model we develop encourages its general use in future investiga- 
tions of nonlinear phenomena. 
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1 Introduction 



Of the many exciting spectacles delivered by Cassini images of complicated 
finescale structure in Saturn's B-ring are some of the most beautiful and in- 
triguing (Porco et ai, 2005). Irregular radial features on large scales, of some 
10 — 100 km, have been associated with the B-ring since Voyager, but Cassini 
has revealed that irregular structure is characteristic of sub-kilometre scales 
as well. In fact, a recent analysis of the UVIS stellar occultation data has 
uncovered quasi-periodic variations on lengths as short as 150 m (Colwell et 
ai, 2007). It is has been suggested that these phenomena correspond to the 
saturated state of the viscous overstability, as examined by Schmit and Tschar- 
nuter (1995, 1999). The viscous overstability is a linear instability that occurs 
when the viscous stresses vary with density in such a way that energy from the 
Keplerian shear can be directed into growing oscillations. The most vigorously 
growing modes possess lengths of some 100 — 200 m and could very well be 
responsible for irregular structure on sub-kilometre scales. On the other hand, 
it is perhaps unlikely that this mechanism generates the much longer 100 km 
features first reported by Voyager, as had been originally hoped. These instead 
may be the responsibihty of ballistic transport (Durisen, 1995) or electromag- 
netic effects (Goertz and MorfiU, 1988). Lending weight to the overstability 
hypothesis are Cassini ISS observations which report that finescale structure 
favours regions in which the optical depth is greater (Figs 5A and 5B, Porco 
et ai, 2005); A^-body simulations show the overstability is similarly sensitive 
to optical depth (Salo et al, 2001). These new observations, and the theoreti- 
cal challenges they throw up, encourage us to examine in greater detail when 
and why this instability occurs (its linear theory) and the manner in which it 
saturates (its nonlinear theory). 

Researchers have mainly probed the onset of the instability with hydrody- 
namic models, sometimes allied with iV-body simulations (Schmit and Tschar- 
nuter, 1995; Spahn et ai, 2000; Salo et ai, 2001; Schmidt et ai, 2001). In 
addition, Schmit and Tscharnuter (1999) have conducted nonlinear hydrody- 
namical simulations and Schmidt and Salo (2003) a weakly nonlinear analysis, 
both with isothermal models. The former study shows that the overstability 
takes the system into a disordered state in which power is transferred to longer 
and longer wavelengths, a process that only saturates if self-gravity is present. 
Schmit 's and Tscharnuter 's adoption of reflecting boundary conditions, how- 
ever, prohibits the development of the stable nonlinear travelling wave so- 
lutions that Schmidt and Salo (2003) discover. If dense regions of the B-ring 
indeed favour the nonlinear wave solution these probably possess a wavelength 
below the resolution of CassinVs cameras but can be captured by the UVIS 
or RSS instruments: in fact, it is very likely that the signiflcant 150 m struc- 
tures reported by the UVIS correspond to these nonlinear waves. Perhaps we 
can then intrepret the irregular structure observed on slightly longer scales 
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(1 — 10 km) as largescale modulations of the wave trains' amplitudes and/or 
wavenumbers, by analogy with water-waves and flame- fronts (Infcld and Row- 
lands, 1990), or alternatively they may correspond to 'defects' between wave 
trains of differing properties, by analogy with the dynamics of the complex 
Ginzberg-Landau equation (Aranson and Kramer, 2002). Certainly there is 
a lot of interesting work to be done examining the nonlinear aspects of this 
problem. 

In this paper, however, we shall reconsider the onset of the viscous overstability 
in a dense ring. We employ a kinetic theory, and hence extend the dilute ring 
analysis of Latter and Ogilvie (2006). In that study it was argued that hydro- 
dynamical models fail to capture important effects associated with the velocity 
anisotropy and non-Newtonian stress. Both we predict will play a part in the 
dynamics of a dense ring as well, a concern which motivates our adoption of 
the better-suited kinetic approach. Moreover, a kinetic model can be framed 
in terms of parameters (particle size, collisional parameters, etc) that can be 
constrained by observation or experiment. This parameterisation also makes 
a comparison with iV-body simulations unproblematic. In contrast, hydro- 
dynamics must make (sometimes speculative) prescriptions for the transport 
coefficients, because their functional dependence on the observables is unclear. 

The price one pays for using a kinetic model is mathematical complexity. The 
dense gas collision terms are especially troublesome as they usually consist of 
complicated multi-dimensional integrals. Consequently, a number of assump- 
tions are necessary to obtain some mathematical purchase on the problem. 
The model we employ is the one developed by Araki and Trcmainc (1986) 
(herafter referred to as AT86), which assumes that the distribution function is 
a triaxial Gaussian. Such a distribution enabled Araki and Tremaine tor educe 
the collision integrals to four dimensions. In this paper we show that further 
integrations are possible even when the coefficient of restitution is allowed to 
depend on normal impact velocity. In fact, if certain additional approximations 
are made, the remaining integrals can be executed by hand. These simpliflca- 
tions greatly facilitate the computation of the equilibrium states (previously a 
difficult numerical problem) , and also open the way for more advanced analy- 
ses: with the elimination of these mathematical obstacles, a kinetic theory can 
be applied more generally, particularly to nonlinear behaviour that is presently 
beyond the range of direct A^-body simulations. 

In many ways, this work is similar to that of Hameen-Anttila who in a suite 
of papers developed the statistical mechanics of orbital elements using many 
of the same approximations and techniques (see Hameen-Anttila and Salo, 
1993, and references therein). The results derived using these methods should 
coincide with ours, though we find that working in a shearing sheet with a 
phase space of velocity and location elucidates local instabilities (such as the 
viscous overstability) much more clearly. 
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The format of the paper is as follows. The rest of this section introduces the 
viscous overstability and sketches out the basic issues pertinent to modelling 
a dense gas. In Section 2 we develop a dense ring model — some of this 
section draws on previous work by Jenkins and Richman (1985) and AT86. 
Section 3 presents calculations and discussions of dense ring equilibria plus 
a comparison with the A^-body simulations of Wisdom and Tremaine (1988) 
and Salo (1991). In Section 4 we undertake a linear stability analysis of these 
states, thereby obtaining criteria for the onset of viscous overstability in the 
non-self-gravitating and self-gravitating cases. These results we compare with 
Salo et al. (2001). In Section 5 we draw our conclusions and point out future 
work. 

1.1 The viscous overstability 

The viscous overstability is one of two instabilities of viscous origin that have 
been proposed to explain structure in planetary rings. The other, the viscous 
instability, occurs wherever the outward angular momentum flux is a decreas- 
ing function of surface density, a situation that initiates a clumping of matter 
into ringlets (Lin and Bodenheimer, 1981; Ward, 1981; Lukkari, 1981). For 
plausible parameters, the angular momentum flux of dense rings has been 
shown to be an increasing function of surface density (Araki and Tremaine, 
1986; Wisdom and Tremaine, 1988), and so it is unlikely that this instabihty 
occurs in Saturn's dense rings. 

The viscous overstability, on the other hand, emerges from a Hopf bifurca- 
tion and, as the name suggests, can be characterised as an overcompensation 
by the system's restoring forces: the stress oscillation which accompanies the 
epicyclic response in an inertial-acoustic wave will force the system back to 
equilibrium so strongly that it will 'overshoot'. Without self-gravity the longest 
lengthscales are the most susceptible to this runaway process, though they will 
grow slowly because the waves' growth rate is proportional to k^, where k is 
radial wavenumber. For sufficently small wavelengths, pressure extinguishes 
the instability and so there is a preferred intermediate scale upon which the 
viscous overstability grows most vigorously. When self-gravity is added the 
disk is more unstable, particularly on a band of intermediate wavelengths. 

The mechanism of overstability relies on: a) the synchronisation of the viscous 
stress' oscillations with those of density, and b) the viscous stress increas- 
ing sufficiently steeply with surface density. In hydrodynamics only the latter 
consideration is relevant, which furnishes the criterion for overstabihty: 

(3 = {d\nu/d\na) > (3*, 

where u is viscosity, a is surface density, and (3* is a number dependent on the 
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thermal properties of the ring (Schmit and Tscharnuter, 1995). GeneraUy, the 
viscous stress need not oscillate in phase with the inertial-acoustic wave, in 
which event the transfer of energy between the two oscillations can be quite in- 
efficient: usually the potential overstable modes are damped as a consequence. 
This is what happens in dilute particulate rings at small and intermediate col- 
lision frequencies (Latter and Ogilvie, 2006). The viscous stress of such rings, 
being predominately local (or 'translational'), possesses a relaxation time of 
order the dynamical time scale, and so they lag behind the epicyclic variations. 
In contrast, the viscous stress in a dense ring is dominated by the nonlocal 
('coUisional') component which oscillates in phase with the inertial motion. 
Moreover, the effective viscosity profiles computed by AT86's dense gas model 
and Wisdom and Tremaine's (1988) particle simulations display a relatively 
steep increase with optical depth (hence a) , suggesting that the effective /3 of 
the system may be sufficiently large to instigate overstability. 

As a consequence, the linear behaviour of the viscous overstability has been 
thoroughly examined, mainly with hydrodynamics (in the papers mentioned 
earlier). However, it is upon the A^-body study of Salo et al. (2001) that we 
will concentrate. There the coUisional evolution of 100 cm radius particles 
were tracked in a shearing box situated in the B-ring undergoing collisions 
that dissipate energy in accordance with the piecewise power law: 



where £ is the coefficient of restitution and Vn is normal impact speed. The 
parameters p and took the values deduced from the ice-collision experiments 
of Bridges et al. (1984), specifically: p = 0.234 and Vc = 0.0077 cm s^^. These 
simulations show that a dense non-self-gravitating disk is stable for all optical 
depths below about four ^ . This critical value is reduced to about one if self- 
gravity's enhancement of the vertical restoring force is included. If self-gravity 
is fully modelled the emergence of the overstability is somewhat retarded by 
the vigorous percolation of non-axisymmetric gravitational wakes, because 
these motions modify the effective viscous properties of the system. 

It is problematic to obtain a general criterion for overstability from simulations 
in terms of the various parameters, as this would require a great number of in- 
dividual runs. Moreover, the fiuid dynamical criterion cannot be neatly related 
to the A^-body results, the main problem being that the quantities involved, 
namely u and are not straightforward functions of the A"-body, or observ- 
able, parameters (particle size, orbital frequency, etc). Some of the hydrody- 
namic quantities can be numerically obtained for a particular equilibrium state 



Throughout this paper 'optical depth' will denote normal geometric optical thick- 
ness. 




(1) 
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using simulations (Salo et ai, 2001; Schmidt et ai. 2001), but the role of the 
various observable parameters is difficult to bring out, as these fluid dynamical 
quantities are tied to the particulars of a specific (and lengthy) computation. 
That said, some of the appurtenant results so obtained (growth/decay rates 
mainly) are in good agreement with the simulations, even if the assumptions 
implicit in a hydrodynamic model possibly lead to errors. 

We seek to use kinetic theory to remedy this problem, by providing a criterion 
for the onset of overstability, specifically an estimate of the critical optical 
depth Tc, as a function of the kinetic parameters of the system, which con- 
nect to the observable data: particle radius a, particle mass m, local orbital 
frequency fl, and the parameters that appear in the elasticity law. 

1.2 Dense gases 

It was not until the late 1980s that it was fully appreciated that nonlocal 
'dense' effects were critical and widespread in Saturn's rings. Perhaps this 
shift was initiated most by the theoretical papers of AT86 and Araki (1988, 
1991), and the simulations of Wisdom and Tremaine (1988) and Salo (1991). 
In fact, the issue had been presaged by Brahic (1977) and examined in some 
detail by Hameen-Anttila (1982) and Shukhman (1984). Essentially, Araki 
and Tremaine and Wisdom and Tremaine revealed that (a) collisional contri- 
butions are crucial to the equilibrium when collisions are sufficiently inelastic 
and (b) the experimentally derived elasticity laws of (1) predict that collisions 
are dissipative enough for this 'cold' regime to be ubiquitous in Saturn's rings, 
given appropriate particle sizes (Marouf et al, 1983; Zebker et al, 1985). In 
addition, both simulations and theoretical models show that the dynamics of 
small particles (for which nonlocal effects are less important) strongly couples 
to the dynamics of the largest particles (for which nonlocal effects are impor- 
tant) both in two-size systems and in polydisperse disks exhibiting power-law 
size distributions akin to Saturn's (Stewart et al., 1984; Lukkari, 1989; Salo, 
1987, 1991, 1992b). 

To capture theoretically this dense regime we turn to the Enskog model which 
revises the Boltzmann theory for moderately dense gases (Chapman and Cowl- 
ing, 1970). The Enskog formalism distinguishes two additional processes not 
captured by Boltzmann kinetics; the first is associated with a large filling 
factor, where filling factor (FF) denotes the proportion of space occupied by 
particles and is defined, for spheres, through 

FF = |7^r^a^ 

where n is volumetric number density and a is particle radius. The second 
process is associated with the collisional transfer of particle properties such as 
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momentum and energy. We describe these in turn. 



When particles take up a significant portion of space (i.e. when FF is not 
small) the volume in which they may move is reduced, and possible collid- 
ers may be screened by other particles (Chapman and Cowhng, 1970). This 
means that the statistics of two impacting particles must include the influ- 
ence of their neighbours and, as a consequence, the evaluation of the collision 
frequency must take space correlations into account. Overall, this leads to an 
enhancement of the collision frequency ujc which, in the Enskog theory, is ap- 
proximately quantified by a factor y(FF) (the 'Enskog factor'). This cannot 
be calculated within the bounds of kinetic theory but must be determined 
separately, usually from data gathered in molecular dynamics experiments 
(Araki and Tremaine, 1986; Jenkins and Richman, 1984; Poschel and Bril- 
liantov, 2004). Inelastic and spinning systems may differ in this respect from 
their elastic and spinning counterparts because they can develop postcoUi- 
sional correlations in their particles' tangential velocity. Inelastic collisions 
diminish the normal component of the relative velocity but conserve the tan- 
gential component, and as a result subsequent collisions lead to a statistical 
alignment of neighbouring particles' traces (Brilliantov and Poschel, 2004). 
Other than possibly forming vortices, this effect can reduce the local value of 
the collision frequency, and so the Enskog factor of an inelastic ensemble may 
be smaller than that of a corresponding elastic system (see Poschel, Brilliantov 
and Schwager, 2002). That all said, in ring applications the net effect seems 
negligible. 

There are two ways properties may be transferred in a particle gas: their 
free carriage by particles between collisions and their transmission from one 
particle to another during a collision. In a dilute gas, the former so-called 
'local', or translational, mode dominates the latter because particles travel 
relatively long distances between collisions. In a dense gas, the mean free path 
is much reduced and can approach a particle radius, in which case the finite 
size of the particles is sufficiently large for the exchange of properties between 
colliding particles to compete with, or even dominate, translational transport. 
Estimates for the relative magnitude of the two types are easy to derive. 
Consider the transport of momentum across a plane by a gas characterised by 
mass density p and velocity dispersion c, and composed of particles of radius 
a. The magnitude of the momentum flux density due to the free carriage of 
particles is of order pc^. On the other hand, the flux density of momentum 
carried in collisions from the centres of all the particles on one side of the 
surface to all the particles on the other side is of order {pc)aujc. This implies 
the scahng: 



where and p* refer to the coUisional and translational pressure tensors 




(2) 
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respectively, Q is the local orbital frequency, and R is defined by 




(3) 



Elsewhere R is called the Savage and Jeffrey i?-parameter (Savage and Jeffrey, 
1981; Araki, 1991), and it quantifies the ratio of shear motions to velocity 
dispersion. Here 17 is a substitute for shear rate. Note that it is (uc/^)R 
that quantifies the importance of collisional transport in the ring dynamics. 
In contrast, an analogous argument shows that R^ quantifies the importance 
of collisional effects in the ring energetics: because heating and cooling are 
proportional to collision frequency the cuc/^ can be dropped from the estimate. 

The effects of large filling factor and collisional transport usually work in 
tandem, though for very high or very low optical depths there are cases when 
one can exist without the other. This is expressed in the scaling: 



For a substantial discussion on this subject see Araki (1991). In Saturn's rings 
R is probably of order unity (Salo, 1992b), which means that filling factor 
effects are as important as the optical depth is large. In low optical depth 
regions, such as the C and D-rings and the Cassini division, filling factor effects 
are hence negligible. In contrast, collisional transport/production effects will 
be important throughout the rings on account of their low velocity dispersion 



Kinetic models and A^-body simulations, in order to make any progress, must 
introduce a number of simplifying assumptions. Typically, ring particles are 
taken to be identical, perfectly smooth, hard, indestructible spheres which 
dissipate energy in coUisions according to a simple 'elasticity law' such as (1) 
(the 'billiard ball' model). Though these presumptions facilitate calculations, 
it is unclear how appropriate some of these are to the real rings of Saturn. 
Neglecting spin and size-distribution should give qualitatively correct results 
for those particles which dominate the dynamics (near 100 cm size). Also it 
is plausible that the timescale of the size dynamics (erosion and accretion, 
etc) may be much larger than the dynamical timescale (il"^), meaning these 
processes do not interfere with the instabilities we study. In any case, the mod- 
elling of spin is hampered by the fact that the frictional properties of the ice 
particles are poorly constrained. We admit that 'smoothing' over the compli- 
cated details of an actual collision with a simple e law is a bold approximation, 
particularly as there is considerable uncertainty about the relevant physics for 
real ring particles. A number of terrestrial experiments, however, have been 
conducted with solid ice spheres (Bridges et al, 1984; Hatzes et ai, 1988; 
McDonald et ai, 1989; Hatzes et ai, 1991; Supulver et ai, 1995; Supulver et 
al, 1997). These show that the coefficient of restitution varies considerably as 



r - FF/R. 



(4) 



(i?~ 1). 



8 



the physical condition of the contact surface is frosty, sintered, or subhmated. 
Depending on the surface condition, various processes ensue, such as erosion, 
sticking, mass transfer, and regohth compactification, all of which influence e. 
But it has also been suggested that a particle may resemble more of a rubble 
pile (a 'dynamic ephemeral body') held together by mutual gravitation (Wei- 
denschiUing et a/., 1984) and/or adhesion (Albers and Spahn, 2006), in which 
case the collisional dynamics will be very different — and the assumption of 
'hardness', and a neat e law itself, perhaps unsupportablc. Unfortunately, we 
lack detailed information about the physical state of a ring particle, and given 
this uncertainty assume that laws such as (1) capture the correct qualitative 
behaviour, though the parameters which appear in them may be subject to 
much variability. 



2 The kinetic model 

2.1 Governing equations 

Consider a gas of indestructible identical nonspinning inelastic spheres of mass 
m and radius a, and with phase space distribution /(x, v, t). The number of 

particles located in the volume rfx centred at x with velocities in the range dw 
centred at v at time t is defined as /(x, v, t)(ix(iv. Number density n, bulk 
velocity u, and the pressure tensor are defined, consequently, through 



The phase space distribution satisfies a kinetic equation which is distinguished 
by a collision operator, {df /dt)c, whose precise form reflects the collisional 
microphysics. This equation is presented after we describe the geometry of 
the problem. 

The particle gas inhabits the shearing sheet (see Goldreich and Lynden-Bell, 
1965). This is a convenient representation of a differentially rotating disk in 
which a small patch, centred on a point moving on a circular orbit at r = Tq, 
is represented as a sheet in uniform rotation, il(ro) e^, and subject to a linear 
shear flow, uq = —2AQxey. The local rectilinear coordinates x and y point in 
the radial and azimuthal directions respectively and Aq = —^rQ^dQ/dr)^. A 
Keplerian disk requires Aq = ^Qq, where flo = fl{ro). This approximation is 
valid if the radial wavenumber of typical variations k satisfies |A;ro| ^ 1, which 




(5) 



(6) 
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is certainly the case for the phenomenon of interest. Prom now the subscript 
'0' will be dropped. 

In terms of the pecuhar velocity, w = v — u, the kinetic equation for particles 
in the shearing sheet is 



df , ,df 



^ - (f ) ' 



where the force per unit mass is 

p ^ d{(t)p + 0d) 



The appropriate centrifugal-gravitational potential of the planet is denoted 
by 0p and the disk's gravitational potential by The tensor eijk is the 
alternating tensor and the angular velocity is $7 = Qe^. By multiplying (8) 
by 1, Wi and WiWj and then integrating over all w we derive the continuity 
equation, 

dfU + dk{nuk) = 0, (9) 

the equation of motion, 

n {dtUi + UkdkUi) = -2neijkfljUk - ndi{(j)p + (po) - djPij + rrii, (10) 
and the pressure tensor equation, 

dtPij + UkdkPij = -PikdkUj - PjkdkUi - PijdkUk - 2eiki^kPij 

- '^^jki^kPii - dkPijk + q-ij, (11) 

where Pijk is the third-order moment, 

Pijk = j WiWjWkfdw, (12) 
the coUisional change in the first moment is 



rrii 



and the coUisional change in the second moment is 



Qij = I WiWj i^] dm. (14) 



c 



For notational brevity we have set di — d/dxi. If nonlocal transport effects 
are neglected then mj is zero on account of momentum conservation. Jenkins 
and Richman (1985) show how to rewrite this term as a flux, so that 
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where pl^ is the coUisional pressure tensor. We will work with this form for 
the rest of this paper. 



The moment equations are vertically integrated, a task that is easily performed 
once we assume that Ux, Uy, and Pij are independent of z. Thus the disk is 
assumed to be vertically isothermal. In addition, vertical averaging permits 
us to close the hierarchy of the momemt equations. In the pressure tensor 
equation Eq. (11), this procedure removes the dominant z derivatives of the 
third-order moment, dzPijk, leaving only its horizontal derivatives, which scale 
relative to the other terms like Hk, where H is disk semithickness and k is the 
wavenumber of a typical horizontal variation. We expect these variations to 
be much longer than the disk thickness and so treat these flux terms as neg- 
ligible. In summary, the closed governing equations wc derive validly describe 
behaviour on scales satisfying H <^ X <^ ro, where A is wavelength. 

We define optical depth as r = na'^N, where N is surface number density. We 
also introduce the velocity dispersion tensor: 

Wij =Pij/n. 

This then allows the vertically-averaged equations to be recast as 

DfT = -rdaUa, (15) 

Dtu^ = -d^i^p + ^d) - 2e^,fsnuf, - - dfs Ww^p + r P^^^/iVl , (16) 
DtWij = -WjkdkUi - WikdkUj - 2ei^k^ Wkj - 2ej^k^ + Qa/N, (17) 

where Dt = dt + u^d^, Greek indices run only from x to y and upper case 
denotes vertical integration except for $p and which are vertically av- 
eraged (if needed). We approximate the central body as perfectly spherical, 
which accounts for the former potential, as then $p = 3Q'^x'^/2 + 0{x/rQ). 
This is acceptable because the precessional effects associated with planetary 
oblateness are unimportant for the instabihties we study. The latter potential 
must be obtained from Poisson's equation: 

V^0D = ^nmGn, (18) 

where G is the gravitation constant 

Finally, we are left with the vertical component of the equation of motion from 
which it is possible to extract an equation for the mean vertical displacement 
of the disk Z (see Shu and Stewart, 1985). In this paper we do not investigate 
vertical warping of the disk. We assume vertical hydrostatic equilibrium and 
symmetry about the plane z — 0; thus Z — P^z — Pyz — 0. The vertical 
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balance furnishes an equation for the vertical stratification of density 



This integro-diffcrcntial equation for n can be manipulated into a second order 
ODE, which must be solved concurrently with Eqs (15)-(17). Alternatively, 
an ansatz for the form of n can be employed, such as a Gaussian (AT86) or 
a polytrope (Borderies, Goldreich and Tremaine, 1985; Mosquiera, 1996), in 
which case (19) becomes an algebraic equation for either the scale height or 
semithickness of the disk, H. 

The assumption of vertical hydrostatic equilibrium, though convenient, may 
be inappropriate when dealing with many dynamical phenomena, including 
overstable oscillations and density waves. This is because the timescale upon 
which vertical equilibrium is established should be of the same order as the 
dynamical time. It is possible to derive a dynamical equation for H, which 
could capture this effect. But for simplicity we persist with the hydrostatic 
model in this paper. 



2.2 Construction of the collision integrals 

In the dilute disk analysis of Latter and Ogilvie (2006) the microphysical 

details of the coUisional interactions could be largely sidestepped by adopt- 
ing a BGK model (Shu and Stewart, 1985; Bhatnagar, Gross and Krook, 
1954). The collision terms in this case were simple, analytic, and local in 
space. Regrettably, no comparable approximation exists for the nonlocal coUi- 
sional interactions of a dense gas. Hence we wade into the mathematics of the 
Boltzmann-Enskog collision theory and build the terms qij and p^j from first 
principles. In their construction it will be assumed that particles only interact 
with each other through collisions and not also through gravity. This omission 
eliminates the possibility of gravitational focusing and scattering. 



2.2.1 The kinematics of binary collisions 

In order to ascertain the consequences of particle collisions collectively we 
must establish the basic physics of a single collision first. Gonsider two col- 
liding particles of mass m and radius a travelling with velocities Vi and V2 
immediately before the collision and v'^ and V2 immediately afterwards. The 
relative velocities pre and postcollision are g = Vi — V2 and g' = v'^ — V2, 
and their locations at the moment of collision are xi and X2. The collision is 
assumed to be inelastic and so the component of g normal to the impact is 




(19) 
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reduced by the factor £, the coefficient of restitution. Hence we write 

g' • k = -£ g • k 

where k is the unit vector directed from the centre of particle 1 to that of 
particle 2 so that 

X2 — xi = 2a k. 

Obviously, g ■ k > (the impact velocity cannot be negative). Also, it is 
assumed that £ is a function of normal impact velocity, v„ = g ■ k. 

This prescription is sufficient to set the postcoUisional velocities: 

v'i = vi-j, V^ = V2+j, (20) 

where 

j = ^(l + £)(g-k)k, 
and mj is the momentum transferred from particle 1 to particle 2. 

The velocity of a particle may be expressed as the sum of the local mean 
velocity u(x, t) and the peculiar velocity of the particle, w = v — u. The 
precoUisional relative velocity is therefore 

g = Wi - W2 + h, 

where h denotes the difference in mean velocity between the locations of the 
particles: thus h = u(xi, t) — u(x2, t). Since the coUision is instantaneous the 
mean velocity does not change in the collision, and because the particles are 
assumed perfectly smooth, the tangential component of g is conserved. 

The prescription above determines the total change in the peculiar velocity 
dyadic in a binary collision, 

A(ww) = w'^w'^ + W2W2 — wiWi — W2W2 

= ^(1 + £)(g • k)[(h - g)k + k(h - g) + (1 + £)(g . k)kk]. (21) 

Hence the total change in peculiar kinetic energy is 

^ mA{w^) = m(l - £^)(g ■ k)^ + m j • h. (22) 

The first term in (22) represents the dissipative loss of energy (zero when 
e = 1). The second term represents the gain in energy from the transfer of 
momentum across the mean shear flow. If the particle size is sufficiently small 
so that |h| -C |wi — W2I (which corresponds loosely to -C 1, see Eq. (3)), 
then this contribution can be neglected. 
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2.2.2 The statistics of binary collisions 

We now determine the frequency of a particular collision specified by vi, V2, 
Xi, and X2. This requires the introduction of the pair distribution /^^^ which 
is defined so that 

/^^^(Vi, Xi, V2, X2, t)(iVi(iV2(iXi(iX2 

is the probable number of particle pairs located in the small volumes rfxi 
and (ix2 centred at Xi and X2, with velocities in the volumes rfvi and dv2 
centred at vi and V2 in velocity space, at a time t. For a dilute gas, in which 
particles are well spaced and travel relatively long distances between coUisions, 
we may appeal to the assumption of molecular chaos (Stosszahlansatz) in 
order to simplify the distribution. This assumption states that there exists 
no correlation between the velocities and positions of any two particles before 
their collision. Consequently, we write /*^^)(vi, xi, V2, X2, t) as the product of 
two single particle distributions, /(xi,Vi,i) and /(x2,V2,t). But in a dense 
gas the probable positions of two colliding particles will be influenced by the 
presence of their neighbours. As mentioned earlier, we may approximate this 
by including the Enskog factor Y, which we take as a function of filling factor, 
and hence number density, n(x). This last assumption limits Y to express 
position correlations only and is justified for elastic gases in which number 
density is nearly uniform. As we only explore homogeneous equilibrium or 
small deviations from it, this aspect of the assumption is acceptable for our 
purposes (AT86). More generally, velocity or higher order correlations may 
need to be taken into account. In summary, we write for two identical particles 
of radius a in contact 

/(2)(vi, xi, V2, xi + 2ak, t) = r [n(xi + ak)]/(vi, xi, i)/(v2, xi + 2ak, t). (23) 

In the limit of vanishing filling factor the Enskog factor approaches unity, and 
we recover the Stosszahlansatz. In the opposite limit, when particles are so 
densely packed that random movement is impossible, Y diverges. This occurs 
when FF = 0.74 for a collection of identical spheres arranged in a face-centred 
cubic lattice (AT86). 

The Enskog factor cannot be calculated from within the compass of kinetic 
theory, but a number of expressions for it have been derived when the kinetic 
gas is elastic, non-shearing, and possessing a small filling factor (Ree and 
Hoover, 1967; Carnahan and Starhng, 1969; Devore, 1984). This paper employs 
for Y the convenient interpolation of AT86 to the molecular dynamics data of 
Alder and Wainwright (Ree and Hoover, 1967). These details appear later. 

We now present the collision frequency. Consider two particles with velocities 
vi and V2 with the first located at x. For a collision to occur characterised by 
a line of centres k in the solid angle (ik, the second particle must be within a 
skew cyhnder of volume 4a^(g • k.)dkdt. Therefore the frequency of a collision 
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between two particles with velocities in the volumes dvi and dv2, centred at 
Vi and V2, with their line of centres within the solid angle dh centred on k, 
and with particle 1 located in the small volume dx centred on x, is 

4 /(2)(vi, X, V2, X + 2ak, i)a^(g • k)dk dvidvadx. (24) 

Into this expression we may substitute Eq. (23) when we are ready. 

2.2.3 The coUisional production integral 

The coUisional rate of change per unit volume of some particle property ip is 
the integral (taken over all possible binary collisions) of the change in in a 
particular collision, multiplied by the probable frequency of such a collision. 
A factor of one half should also be included, otherwise we will count every 
collision twice (once for particle 1 and once for particle 2). We denote the 
resulting quantity by C{ip) and express it as 

C{^)=2 1 i/(g.k)A(V')/('Hvi,x,V2,x + 2ak,t)a2(g-k)dkdvidv2, (25) 

where the velocity integrations are over all space, the k integrations are over 
the surface of the unit sphere, and H is the Heaviside step function (Jenkins 
and Richman, 1985). The step function ensures that the integration only in- 
cludes real coUisions, i.e. ones that satisfy g • k > 0. Note importantly that 
C{ijj) is nonlocal in space, because production of at x necessarily depends 
on the values that certain fields take at small but finite distances away. To 
determine the coUisional production of second moment, ip = ww, we use the 
expression (21) for A(ww). We subsequently define q = C(ww) (cf. Eq. (14)). 

If the system is dilute we set FF ^ 1 which forces "K ~ 1. Also we do not 
distinguish between the centres of the two particles, i.e. we assume the mean 
flow and the distribution function do not vary significantly on length scales 
of order a. Thus h and /(v, x + 2ak, t) ^ /(v, x, t). If we set i/j — ww, 
these simplifications reduce expression (25) to that of Eq. (24) in Goldreich 
and Tremaine (1978). 

2.2.4 The coUisional flux of momentum 

Consider the momentum transferred in a collision from particle 1 to particle 
2 when they both straddle a surface of area dS with normal n. The normal 
component of momentum transferred across the surface from particle 1 to 
particle 2 is then mj ■ n, and the condition that the particles straddle the 
surface in the right sense is k • n > 0. 

Consider all the momentum transferred across the surface by such collisions. 
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Particle 2 must lie just above the surface and sit in a skew cylinder of volume 
a(k-n)rfS'. But the particle striking it must lie in the cylinder 4a^(g-k)(ikdt (for 
collisions within solid angle and during time dt). Therefore, the average 
collision rate for particles with velocity ranges dvi and dv2, straddling the 
surface dS, oriented in the range dk with particle 1 at x is 

4/(2) (vi, X, V2, X + 2ak, t)a^{k • n)(g • k)dkdvidv2dS, (26) 

provided g • k > and k • n > 0. The total momentum flux density resulting 
from collisions we denote by and deflne through 

n-p' = Ama^ J H{g- k)/^^) (vi, x, V2, x + 2ak, (k • n) (g • k) j dk dvi dv2, 

thus 

p" = 2ma^y i/(g-k)/(2)(vi,x,V2,x+2ak,t)(g-k)2(l+£)kkdkdvidv2. (27) 

If we assume that the mean flow u varies slowly across a particle diameter, 
then we may write 

h = u(x, t) — u(x + 2ak, t) ^ —2a k ■ Vu, 

in which case we have the integral equivalent of Eq. (22): 

^mtr(q) = -L> - : Vu (28) 

where the rate of dissipation is 

D^^a'^mJ H{g ■ k)/^') (vi, x, V2, x + 2ak, t){l - e^) (g • kfdkdvidv2. (29) 

This relates the coUisional production of energy to the collisional momentum 
flux, and shows that the formalism is consistent; the quantity p'^ : Vu is just 
the rate at which the collisional momentum flux extracts energy from the 
mean flow. 



2.2.5 The peculiar velocity distribution function 

To make any progress in analysing the governing equations the integral ex- 
pressions must be simphfled. To that end we propose a speciflc model for 
the distribution function. First we introduce the peculiar velocity distribution 
function /^(w, x, i) deflned so that 

/(v, X, t) = n(x, t)/w(v - u(x, t), X, t). 
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Following Goldrcich and Tremaine (1978) and AT86, atriaxial Gaussian model 
for /w(w,x, t) is adopted: 



U = [(27r)'det(W)]-i/2gxp 



-iW-^ : WW 



(30) 



where the dependence of /w and W on x and t is understood. This form 
satisfies the required moment properties of Eqs (5)- (7). 

To facilitate the manipulations of the next section we now assume W does 
not vary appreciably on lengthscales of order a. This is acceptable if we only 
wish to analyse variations on long radial scales. In fact, if a ^ H then the 
assumption is consistent with the closure scheme in Section 2.1. Note that we 
have already assumed vertical isothermality, and so W's variation with z is 
neglected. 



2.2.6 Transformation of the integrals 

First we transform from the variables (vi,V2) to the relative and centre of 
mass velocities, (g, Vc) , the latter defined by 



Vc = -(Wi + W2). 

The Jacobian determinant of the transformation is unity and the only appear- 
ance of Vc in the integrands is in 



/w(wi) /w(w2) = [(27r)Met(W)J exp 
This allows us to integrate over Vc- 



-W-i:(vcVe + |(g-h)(g-h)) 



/ 



exp 



7r3det(W)] 



1/2 



What remains is 

q^2a'J H{g • k)(g ■ k) F(x, k)A(ww)/g(g - h) dkdg, 
= 2a'm J H{g • k)(g • k)^ F(x, k) (1 + e)f^{g - h) kkdkdg, 

where we have introduced, for notational brevity, the function 

y(x, k) = y [n(x + ak)] n(x) n(x + 2ak), 
and a distribution function for g, 

/g(g) = [(47r)=^det(W)]"'^%xp (-^W"^ : gg) . 



(31) 
(32) 

(33) 
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The latter is normalised such that / fgdg = 1. Note that /g(g — h) is an off- 
centred triaxial Gaussian, and the offset h which derives from the background 
mean flow, will depend on k. 

We now tackle the g integrals. We express g space in Cartesian coordinates 
{X,Y,Z), orienting them so ez = k. Thence we can remove the step func- 
tion and replace the integration domain with the half space, {{gx, Qy, gz) '■ 
gx, gr ^ gz G M+}. This choice also simplifles the functional dependence 
of the coefficient of restitution to £ = e{gz)- Consequently, the gx and gy 
integrations are straightforward, and the gz integrations go through if we in- 
troduce certain 'averages' of e. The biggest problem is returning the result to 
coordinate- free form, which requires some algebraic insight. But once this is 
done we have 

q = 27r-VVyyC'/'x 

{ C( (1 + )3 Fi K - [ (1 + {s)s) F, - e(l + {8)2) F2 ] G } dk, (34) 
p-^ = 27r-i/2^3 JyC{1 + (6)2) F2 K dk, (35) 

where the dyadics K and G are defined through 

K = kk, G = (Wk)k + k(Wk), (36) 

the F functions are 

Fi = 2(1 + e)e-^'' + TT^/^e (3 + 2eKl + erf 0, (37) 
F, = 2ie-^' + 71^2(1 + 2e){l + erf 0, (38) 

where erf is the error function; the two variables C and ^ are defined through 
and the average by 

for some general function t\){x). If 7/; is a constant then {i))p = ip. Expressions 
(34) and (35) are analogous to Araki and Tremaine's (1986) but are more 
general as they include the dependence of e on impact velocity. They are sim- 
pler, requiring the evaluation of two-dimensional rather than four-dimensional 
integrals, and also benefit from being coordinate free. 

The quantities that appear in Eq. (39) have straightforward interpretations: ( 
quantifies the magnitude of the squared impact velocity for a given orientation 
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k, while ^ measures the ratio of the shear velocity across a particle radius to 
the velocity dispersion, and is hence oc R. This means dilute expressions may 
be recovered by letting ^ = and Y = n(x)^. 

In order to evaluate the remaining integrals, expressions for e and h must 
be supphed. The experimentally derived, piecewise power law of (1) is the 
obvious choice for e, while we treat the velocity gradient as locally uniform, 
as before, i.e. h = —2a k ■ Vu. Therefore: 



where e is the rate of strain tensor. 

In addition, we require the form of n(x) and y[n(x)]. The Enskog factor 
may take one of the formulae mentioned, but the number density should be 
determined from the vertical momentum balance. Alternatively, an explicit 
form could be supplied. We discuss these issues a little later. 

The integral terms as they stand, though still involved can be computed nu- 
merically without much difficulty. Nonetheless, additional assumptions yield 
further useful simplifications. 

First, we may suppose that n does not vary appreciably on scales of order a; 
so, like W, it may be considered to be the same at the position of each particle 
in contact: i.e. Xi = X2 = x. Then 



and can be brought outside the integral. This means that as far as the collision 
frequency is concerned nonlocality is neglected and n is constant; for all other 
purposes (in the height equation, etc) n varies and must be determined from 
(19). If we consider homogeneous equilibria or largescale horizontal variations, 
n will only exhibit appreciable variations in the z direction (the disk thickness 
being of order a) and so, essentially, it is this vertical variation that Y omits. 
We hence refer to the adoption of (41) as 'vertical locality', as opposed to the 
'vertical nonlocality' of the full model (Eq. (33)). Note that nonlocality in the 
shear, introduced by h, remains and is crucial. 

Second, we may suppose that for the purposes of the colhsional kinemat- 
ics the coefficient of restitution is an averaged quantity dependent on only 
macroscopic variables, principally c, as was done in Latter and Ogilvie (2006). 
Hence, e is constant as far as the integrations are concerned. 

Third, we may approximate the integrands by polynomial interpolation or se- 
ries expansion. These approaches are valid if dense effects and/or anisotropy is 



a 




(W : kk)V2 



Y = y[n(x)] n(x)2 



(41) 
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small, and permit the remaining integrations to go through. We shall examine 
the last two assumptions in some detail. 



2.2.1 Averaging the coefficient of restitution 

As we have seen, the expressions for q and p'^ involve various averages of 
e which are described in Eq. (40). These can only be analytically evaluated 
for relatively simple functions, and then the result is rather complicated. For 
example, if we set the coefficient to be a straight power law, £ = (g • k/vc)""^, 
we get 



(->-^r(2-ip)r(i-irt^^j .^_jH,_3(_«, (42) 

(e>3 ~ r(| - ip)r(2 - Ip) (^) ' . (ij) (43) 

where H,y is the Hermite function, F is the gamma function (Abramowitz and 
Stegun, 1965), and Fi and F2 are given by Eqs (37) and (38). Unfortunately, 
more complicated functions, including the more accurate piecewise law, do 
not yield a closed form for {e)g, in which case the collision terms require the 
numerical evaluation of three integrals each. 

We aim to employ the piecewise e law and, given that the collision terms are 
sufficiently complicated, it would be convenient if we could eliminate the third 
integral. Moreover, expressions like (42) and (43) are complicated functions of 
the integration variable k and hence computationally more intensive. These 
consideration motivate the technique of 'preaver aging'. Because the elasticity 
laws we possess are in some sense 'the average' of a number of experimental 
collisions anyway, philosophically, at least, this approximation is justified. In 
this subsection we provide details of a possible preaveraging process. 

Let us define A^{ip) as the weighted average of the function ipi^g ■ k) over all 
collisions. The weighting function is (g-k)'''""^, and 7 is a real positive number. 
We define 

, . / H{g ■ k) ^(g ■ k) Y hh (g ■ k)-^ <P\<id^w,<Pw2 

^-rW) - JHig.k)Y hh (g • k)7 d?M^^,<Pw2 ^^^> 

where /i and /2 represent the peculiar velocity distribution function evaluated 
at particle 1 and 2 respectively. Obviously, A^{'^) = -0 if '0 is a constant. We 
now replace each instance of {e)q with Aq{e), i.e. {e)q ^ Aq{e). For a simpler 
but less accurate formalism, we replace every instance of {e)q with A^{e) for 
a fixed 7 for all q. 

The integral (44) can be evaluated by the same techniques as earlier if it 
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is assumed that / is a triaxial Gaussian. However, the result is relatively 
complicated and it so happens there is little gain over the fully consistent 
model of (34) and (35). An acceptable approximation is the adoption in (44) 
of a MaxweUian, / oc exp [—w'^/{2c^)]. However to complete all the integrals 
we must adopt (41) and neglect the nonlocal contributions, i.e. set ^ = 0. This 
is less justifiable but mathematically advantageous and, in fact, yields results 
in good agreement with the full model, as we explore in Section 3.10. 

Implementing these approximations for the piecewise power law function, e — 
min{l, (g • k^c)"^}, yields 

A2(s)^eriX-^X (l + —xA e-""" - X' E^,{X^) (45) 

As{e)^l-(l + X'-- xA e-^^ - 7^ E,, {X') (46) 

where X = Vc/{2c), fig = (p — 1 — s)/2, and E^{x) is the exponential integral 
E function. Though they may seem a little complicated, these expressions are 
especially helpful because they involve only the squared velocity dispersion 
(c^) and none of the integration variables (through ( and ^). Consequently, 
they can be brought out of the integrals entirely. 

2.2.8 Interpolating the integrands 

Because the F functions rapidly increase with ^, their Taylor series (trun- 
cated at a manageable order) are poor approximations unless ^ is very small. 
But dense ring equilibria fall into the regime of \^\ ~ 1. An alternative is to 
interpolate the /'s on some range of ^, which we can estimate from Eq. (39). 

To proceed, vertical locality is assumed, so that Y = F(x)n^(x), and e is 
treated as preaveraged, with the additional approximation, ^2(5) = A^i^e) in 
the coUisional production integral. Enforcing these gives 

q=27r-V2a2f[l + >l3(£)]|[l+A3(£)] J C^/^FiK dk - y C'/^FaG dkj , 
= 27r-^/'a^y[l + A2ie)] J ( F2K dk, 

where F3 = Fi — ^F2. 

We now interpolate the F functions. For small anisotropy we have the upper 
bound 





e 


\C\<-R 









which for equilibrium solutions with R < 1 gives |^| < 9/8. On this range we 
find that Fi can be successfully interpolated by a cubic in ^, and F2 and F3 
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by quadratics: 

ra=0 n=0 n=0 

where the Fjn are pure numbers arising from the interpolation range and the 
interpolation method. 



Next, we express 



W = c^il + a) 



where the squared velocity dispersion is = Tr(W)/3 and a is the anisotropy 
tensor, which is symmetric and traceless. It follows that 

TT = —mn (? a 

if the viscous stress is denoted by tt. Because ||a|| < 1, the integrands may be 
systematically expanded in a through 

e = (e : kk) [ 1 - ia : kk + |(a : kk)^ + ^(a^) 

and 

C = [1 + a : kk] . 

Now we are in a position to evaluate the integrals directly. Expressions for q 
and p*^ may be written in terms of spherical surface integrals of the form 



where S is (the surface of) the unit sphere. Necessarily 

1 , 



'Ill2---l2l 21 -\- 1 

where 

Ihi2—i2i — ^{1112^1314 ■ ■ ■ ^i2i-ii2l)y 

and is the unique, totally symmetric, isotropic tensor of rank 21. Parantheses 
around a set of j indices denote the sum of all permutations of the indices 
divided by j\ 

The closed forms of the integrands to 0{sl) are presented in Appendix A, and, 
though complicated, are algebraic functions. As such, they possess a great 
numerical advantage over previous kinetic approaches, and certainly pave the 
way for dynamical analyses, particularly nonlinear studies, which have been 
the sole preserve of hydrodynamics. 

The main errors that characterise this formalism arise from truncation in a. 
Thus, significantly anisotropic systems may not be well described (for instance. 
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rings suffering very few collisions). In addition, when the velocity dispersion 
is extremely low or the rate of strain large, the interpolation range must be 
extended, and at some point the interpolating polynomials may be of too low 
an order to adequately approximate Fi, F2, and F3. This should only be an 
issue for some nonlinear calculations, however. 



2.2.9 Expanding the integrands 

We have argued that a real ring lies in the regime ^ ~ 1, but it is nevertheless 
interesting to examine 'warm' and dilute rings for which ^ is small. Doing so 
allows us to connect our formalism to other well-known kinetic models, such 
as those of Goldreich and Tremaine (1978), and Jenkins and Richman (1985). 

We expand the integrands of q and p'^ first in powers of ^ and then in a, 
assuming both are small. All the integrals then can be performed much as 
before. In the interest of space we omit these details. 

The closed algebraic expressions for q and p'^ so obtained are analogous to, 
and possibly more systematic than, those computed by Jenkins and Richman 
(1985), who employ Grad's prescription for the single particle distribution 
function (Grad, 1949). They concentrate on the case when the spatial deriva- 
tives of the mean fields are small, and the distribution nearly isotropic. In 
the language of our model, they examine orders up to and including (,^°, a^) 
and (,^^, a'^), and then in the hydrodynamic limit, oOc large. This permits them 
to obtain constitutive expressions for the hydrodynamic transport coefficients 
— shear and bulk viscosities and thermal conductivity. Interestingly we can 
recover precisely the same expressions, despite the difference in our choice of 
distribution function. 

The leading terms, of order are just the coUision integrals for a dilute gas. 
And because of the triaxial Gaussian ansatz adopted they coincide with the 
collision term of the Goldreich and Tremaine (1978) model up to whichever 
order in a we choose to truncate. After vertical averaging we obtain 

Q = ^ T^Nc^il + e) - e)l - i(3 - £)a + 0{s?) 

In contrast, the model of Shu and Stewart (1985) gives 



i(i-)i-(T^, 



and their interpretation of the 'dilute' Hameen-Anttila model (Hameen-Anttila, 
1978, 1981) is 

Q^^TnNc^{l + e) [-|(1 -£)1 - 1(3 -£)a 
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Here upper case indicates vertical averaging, as before. All these models agree 
at (9(a°) but disagree at 0(a). However, the difference at this order between 
the Goldreich and Tremaine and Hameen-Anttila models is only a small dif- 
ference in the constant coefficient (1/6 against 1/5). This accounts for the 
(relatively) good agreement these two models enjoy in equilibrium calcula- 
tions (Stewart et al, 1984). 



3 Equilibria 

Now that the various models have been sketched, we employ them to calculate 
dense ring equilibria. We solve for the vertically averaged pressure tensor and 
the disk height for a given vertical structure of n. First, we compare the 
results predicted by our kinetic model with those produced by the kinetics of 
AT86 and the A'"-body simulations of Wisdom and Tremaine (1988) and Salo 
(1991). Second, we test the appropriateness of the three models introduced 
in the preceding sections, namely, that of 'vertical locality', the 'preaveraged' 
e model, and the interpolation model. Because these formalisms significantly 
simplify the problem it is worthwhile testing the approximations upon which 
they are founded. 

3.1 Governing equations 

We calculate the steady homogeneous state of a disk rotating about a perfectly 
spherical planet and whose self-gravity is negligible in the horizontal balances. 
In the shearing sheet this equates to n = nQ{z) and u = —{3/2)flxey. The 
equilibrium pressure tensor equations are 



where, as earlier, upper case indicates vertical averaging (except in the case 
of Wij and c which we assume independent of z) . The natural nondimension- 
alisation of this system is founded on particle radius and orbital frequency, so 
that 



where asterisked quantities are dimensionless. 

In addition to this set of equations we need expressions for nQ{z), Y[n] and e. 



^ nW^^ + 29Wyy + Q^y/N 
-nW^y + Qyy/N 

Qzz 



0, 

0, 



(48) 
(49) 

(50) 
(51) 



t = n-H\ x = ax\ u = anu*, Wij ^ a^n'^W*., 
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The last two can be supplied by simulations or experiments. The first can be 
determined from the vertical hydrostatic balance (19). We visit each in turn. 



3.2 The Enskog factor 

The Enskog factor is modelled according to Araki and Tremaine's prescription: 
a polynomial fit to the molecular dynamics data of Alder and Wainwright (in 
Ree and Hoover, 1967) ^ : 

5 

y(FF) = Yl AFF^ (52) 

i=0 

where 

Do = 1.0000, Di = 3.5496, D2 = -18.816 

D3 = 165.30, D4 = -407.27, D5 = 406.94 

Though strictly applicable to only non-shearing gases of elastic spheres, this 
prescription should capture the qualitative behaviour of the Enskog factor in 
an inelastic system (see AT86). 



3.3 The coefficient of restitution 

For e we mainly use the piecewise power law of (1) which has been experimen- 
tally derived by Bridges et al. (1984), Hatzes et al. (1988), and Supulver et al. 
(1995) for sufficiently frosted ice particles in Saturnian conditions. The two 
parameters p and Vc we let vary, though recognising that Bridges' data give 
their values as 0.234 and 0.0077 cm s^^, Hatzes' data 0.20 and 0.025 cm s~^, 
and Supulver's 0.19 and 0.029 cm s^^. The preaveragcd e's take the functional 
form of (45) and (46) in the collision terms. For one calculation we employ 
the power law form of e. for which wc use the fully consistent expressions 
(42) and (43). These forms introduce the important dimensionles quantity 
Rc = aQ/vc which controls to some extent the equilibrium velocity dispersion 
of the particles. 



^ It should be noted that the values given by Araki &: Tremaine (1986) contain 
a typographical error. Their coefficient -D3 must read 0.8874 not 0.08874 (Juergen 
Schmidt, private communication). 
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3-4 The height equation and vertical locality 



Instead of solving for n a convenient alternative is to assume it takes a specific 
form (AT86), either a Gaussian or polytrope. In this case Eq. (19) becomes 
an algebraic equation for H. This is then solved with the set (48)-(51). 



In the case of Gaussian vertical structure, we write, 



no (2;) = n exp 



^2 



2H^ 



(53) 



where H is the disk scale height and n a constant. Consequently we can 
determine the optical depth from 




where FFq = FF[n], i.e. the filling factor in the midplane {z — 0). This 
expression gives flesh to the estimate (4) . 

After multiplication by z and vertical integration, the hydrostatic balance (19) 
becomes 

n^H^[l + lV2 Gs FFo] - W,, - Pl/N = 0, (55) 
where N — ^/2n Hn is the surface number density and 

expresses the disk self-gravity. Here Pp and ps are the internal mass den- 
sities of a particle and the planet respectively, r is distance from the cen- 
tre of the planet, and rs is the radius of the planet itself. In Saturn's case, 
Ps = 618kgm~^ and rs — 6.03 x 10^ km; we then have 

= 0.735 xf J'^ ){ -]. (56) 

VlO^km/ \^100kgm-3y ^ ^ 

For particles composed of solid ice, Pp = 900 kg m^'^, the parameter Gg varies 
from about 6 at the inner edge of the B-ring, to 7.7 at the outer edge of the 
B-ring, to approximately 9 at the outer edge of the A-ring. Of course, these 
estimates are rather rough. The internal density of real ring particles may 
depart signficantly from that of solid ice: particles may, for instance, be far 
more 'fluffy' and possess a density less than that quoted (see, for example, 
WeidenschiUing et al, 1984). 
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A polytropic model can be described by 



(2 \ 
1 - (57) 

for a given constant a. Then the optical depth is 

^- 4r(3/2 + «) 7^^° ^^^^ 

and the vertical hydrostatic balance differs httle in form from (55) except for 
the addition of the new parameter a, 

Ql (\{a + |)-i + I g{a) FFq) - W,, - P'JN = 0, (59) 

Vi(a + l)'r(2a + |)- 
This approximation has been used by Borderies, Goldreich and Tremaine 
(1985) and Mosquiera (1996) (with a — 1) and is more accurate than the 
Gaussian model at describing monodisperse disks at moderate and large fill- 
ing factors, as we shall see. A problem, however, is the specification of a 
which may vary with FFq but which cannot be determined from the govern- 
ing equations. That said, our equilibrium solutions show that the equilibrium 
quantities are relatively insensitive to the precise choice of a. 



where 



If we assume vertical locahty, the scale height H can be removed from the 
viscous stress equations (48)-(51) completely. Moreover we can solve for H 
explictly in (55) or (59) once the viscous stress components are known, because 

then Y/N, and thus P^^/N, is independent of H. In the case of the vertical 
Gaussian, the averaged Y is 



>'/'V=J^(g«,FF»VV.+2|. (60, 

In the case of the polytrope, we compute 

Y/N (yd- FFt^ r(a + 3/2)-r(ic. + l) . 

""^'^ IGnaA^,''"' r(a + l)-r(ja + 3/2) '• ^^'^ 



If we do not adopt vertical locality then H appears in all the vertically averaged 
equilibrium equations through the vertically averaged Y (cf. Eqs (34) and 
(35)). The height equation must then be solved concurrently with (48)-(51). 
Unfortunately, the vertical integration can only be accomplished analytically 
for the Gaussian model. Usage of the polytrope model must be accompanied 
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by the assumption of vertical locality, unless we compute the extra integral 
numerically. 

In the Gaussian case, vertical nonlocality engenders the additional factor 
exp [— (a^/i7^)fc^] in (60). This 'weighting function' diminishes the effect of 
those collisions in which the particles inhabit different vertical strata, or, put 
another way, arc oriented 'more vertically' with respect to each other (ex- 
hibiting larger kz)- It follows that the factor takes values between 1 ('perfectly 
horizontal' collisions) and exp[— a^/i7^] ('almost vertical' collisions). When 
the disk is thick (meaning, a/H is small) then this range is very narrow and 
the factor is essentially 1; vertical nonlocality is then unimportant. In the ex- 
tremely thin limit, i.e. a/H ~ 1, the factor can be smaller than e~^, which 
expresses the simple fact that in a near monolayer it would be unusual to find 
a pair of particles inhabiting very different vertical positions. 

Finally we present an expression for the collision frequency. This should scale 
roughly like na^ cYn. After inspecting (24), we settle on 



in which the constant coefficient is chosen so that the expression conforms 
with the dilute ujc calculated by Shu and Stewart (1985). In this limit the 
vertical structure is Gaussian and we obtain a simple linear dependence on 
optical depth: Uc = {8/n)QT. 



3.5 Numerical technique 



The equilibrium solution requires the roots of the algebraic system (48)-(51); 
if vertical nonlocality is included we must add the the height equation (55) or 
(59) to this set. As these equations involve integral functions of Wij a number 
of numerical techniques arc needed. The solutions are obtained by using the 
multi-dimensional Newton's method, while the spherical surface integrations 
are broken into azimuthal angle, which are accomplished using the trapezoidal 
rule, and meridional angle, which are accomplished using Gaussian quadrature. 
Mathematica was the programming language used, and the special functions 
that appear in the integrals are evaluated using the Mathematica routines. If e 
is specified, as either a constant or a function of g • k, the equilibrium solution 
Wij (and H) depends only on the parameter FFq. 




(62) 



28 



3.6 Comparison with the model of Araki and Tremaine (1986) 



In order to check our model we compared its predictions with those of AT86, 
who adopt a constant e, Gaussian vertical structure and vertical nonlocality. 
If we do the same, the two models should produce identical results because 
they share exactly the same assumptions. 

We computed equilibrium properties with our formalism and checked it against 
those produced by the Araki and Tremaine model (Schmidt, private commu- 
nication). The parameters used were e = and Gg = 7.4, on one hand, and 
e = 0.5 and = on the other. Agreement to four significant figures was ob- 
tained and this seemed only limited to the accuracy chosen for the numerical 
integrations. These results prove that the formalism we develop in this paper 
is equivalent to that of AT86. 




Fig. 1. In these figures we show the effect on certain equiUbrium properties when 
a factor 2 is introduced into the colUsional transport term, P^^. Additionally, we 
plot points from Fig. 2 in AT86 in order to demonstrate that their computational 
routine suffers the factor 2 error. In subplot (a), we present the orientation angle of 
the velocity ellipsoid 5 and the vertical velocity dispersion C3. The parameters taken 
are e = and Gg = 7.4. In subplot (b), we present optical depth r. The parameters 
are e = 0.5 and Gs = 0. Generally the discrepancy is small, except in r and H when 
self-gravity is weak. 

It should be mentioned that the kinetic equilibria in Araki and Tremaine 
(1986), and also Wisdom and Tremaine (1988), were computed by a routine 
which introduces an erroneous factor of two in the calculation of P^^ (given 
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by their Eq. (142)). As a consequence their computations overestimate the 
amount of colhsional momentum transport in the vertical direction. In the 
self-gravitating case Gs = 7.4, which appears in AT86, the discrepancy is 
small but noticeable. We illustrate the discrepancy in Fig. la and also plot 
some points from Araki and Tremaine's Fig. 2. As is clear, the latter sit per- 
fectly on the calculation that doubled P^^. In contrast, when self-gravity is 
omitted, the overestimated collisional flux has greater influence, but only on 
those quantities which depend closely on the vertical problem, namely optical 
depth r and scale height H. This is demonstrated in Fig. lb where we plot 
T versus FFq for the case e — 0.5 and Gg — 0. For large FFq the discrepancy 
can be as large as 35%. 



3.7 Comparison with the simulations of Widsom and Tremaine (1988) 

The discrepancy just mentioned means that the agreement between the sim- 
ulations of Wisdom and Tremaine (1988) and the AT86 kinetics with respect 
to T no longer holds, cf. their Fig. 15. We stress, though, that the good agree- 
ment witnessed by the 'horizontal' equilibrium properties, (as functions of 
FFq) remain. 

The disagreement in r we blame on the assumption of the Gaussian vertical 
profile, which is appropriate at low FFq but less so at higher filling factors, 
when the disk establishes a vertical structure akin to a polytropc (see Fig's 
3, 4, 8, and 9 in Wisdom and Tremaine, 1988). These two regimes and the 
transition between them can only be captured by a kinetic model that properly 
solves the vertical hydrostatic equilibrium, which is precisely what is done in 
Araki (1991). The agreement with Wisdom and Tremaine's simulations in this 
work is excellent, though the kinetic theory is computationally expensive and 
suffers numerical difficulties. 

As an alternative to the approach of Araki (1991), we apply a model in which 
polytropic vertical structure is assumed from the outset. Such a model should 
approximate well the middle to larger filling factor regimes and tolerably ap- 
proximate the low filling factor regime, thus saving us the arduous task of 
computing the vertical hydrostatic balance exactly. 

Subsequently, we computed equilibria using the polytropic formalism embod- 
ied in (57)-(59), and (61) using a fixed a and e — 0.5 and Gg — 0. Much better 
agreement was obtained with the Wisdom and Tremaine simulations at mid- 
dling and larger filling factors and reasonable agreement at lower filling factors 
(cf. Fig. 2). Indeed, Fig. 2a nicely illustrates the transition between the two 
regimes — Gaussian at lower FFq and polytropic at middle to large FFq. The 
transition occurs around FFq = 0.275 or about t — 1. Additionally, we plot 
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Fig. 2. (a) Here we plot optical depth function of central filling factor FFq, 

computed by various methods. The dashed and dotted lines correspond to kinetic 
models with polytropic vertical structure ('PT') and a = 0.4 and a = 1 respectively. 
The solid line corresponds to a kinetic model with Gaussian vertical structure. The 
diamond points (WT88) correspond to the simulations of Wisdom and Tremaine 
(1988), read from their Fig. 15. (b) The second panel presents velocity dispersion 
C, as a function of central filling factor, computed by the three kinetic models. The 
coefficient of restitution is set to 0.5 and = in all cases. 



the velocity dispersion versus FFq to show, in contrast, that the horizontal 
equilibrium properties are little altered by the choice of vertical model. 



Fig. 2 



The range of a that we employed was informed by simple inspection of the 
vertical profiles of Wisdom and Tremaine (1988). These indicated that a takes 
values between 1 to about 0.4 (and perhaps lower). It is clear, though, that 
changing a does not significantly alter the equilibrium properties, certainly 
not for Wij, as can be seen in Fig. 2b. Thus it appears the polytropic model 
constrains the equilibrium behaviour to a separate regime effectively indepen- 
dent of the specific choice of a: this means we can't just produce any result 
we like by tweaking a. However, the best agreement is obtained with a = 0.9 
and we employ that value throughout the rest of the paper. 
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3.8 Comparison with Salo (1991) 



We next compare our equilibrium results with those of Salo's (1991) simu- 
lations for which the coefficient of restitution varies as a function of impact 
velocity. This supplies us with results with which to test the kinetic model's 
treatment of nonconstant e. Moreover, this set of data is particularly useful 
as it includes a number of runs with different fc, one of the key parameters 
in our model. The disk, otherwise, is non-self-gravitating, and particles are 
assumed to possess a radius of 100 cm and orbit at the location in the B-ring 
where aVt = 0.02 cm s~^. Collisions dissipate energy according to the piece- 
wise power law, Eq. (1), with p = 0.234 and Vc a multiple of Vb = 0.01 cm s^^ 
By varying Vc the parameter Rg can be controlled, which roughly governs the 
velocity dispersion and the 'denseness' or 'diluteness' of the gas, at least for 
intermediate collision frequencies. 

We will actually compare our results with more recent equilibrium simulations 
which reproduce the runs of Salo (1991) but with a five-fold increase in particle 
number and duration (Salo, private communication). The results are much the 
same but with reduced scatter. 

The kinetic model adopted in this section employs 

• a preaveraged e, the piecewise power law, using (45) and (46) 

• polytropic vertical structure with o; = 0.9 

• vertical locality. 

When nondimensionalised the model depends on p, = aQ/vc, Gg, and 
FFq. In this section = and p = 0.234 as mentioned, and FFq and Rc = 
2(ff)/fc) will vary. Once supplied with these inputs our algorithm produces the 
nondimensionalised Wij and then computes a number of related properties: the 
viscosities, 

•in/2' 3Q/2 ' 

the disk semithickness i/, optical depth r, and preaveraged coefficient of resti- 
tution, A^{e). Throughout this section the quantity Ay,{e) will be referred to 
as simply e. 

3.8.1 Results 

In Fig. 3 we plot the equilibrium velocity dispersion, the disk scale height, the 
averaged e, and the central filling factor against r for various Vc- In Fig. 3a we 
also include the results of Salo's simulations. 

As is clear from Fig. 3a the kinetics and the simulations are consistent, and 
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match particularly well when filling factor is high. In general, the kinetic 
model slightly overestimates the velocity dispersion, particularly at low opti- 
cal depths. The preaveraging assumption probably contributes most heavily to 
this discrepancy. Errors arising from the asumptions of vertical locality should 
be secondary because the disk thickness H is relatively large (cf. Fig. 3c). But 
the choice of polytropic vertical structure may play some part, even though 
Fig. 2b suggests that a Gaussian model would also overestimate C. For large 
ratios Vc/vi^ — when the velocity dispersion is large and FFq low (see Fig. 3d) 
— the vertical disk structure will probably resemble a Gaussian not a poly- 
trope. 

These considerations aside, the general trend is quite clear: as Vc increases so 
does the equilibrium velocity dispersion. This admits a relatively straightfor- 
ward interpretation. Because coUisions must dissipate energy to offset viscous 
heating, C must be at least 0{vc): if C were much less than Vc very little en- 
ergy at all would be lost in collisions and certainly not enough to maintain an 
energy balance. It must follow that if Vc is increased so must the equilibrium 
C. 

The magnitude of the velocity dispersion has a direct bearing on the thickness 
of the disk as Fig. 3b can attest. Here we can see how the particles' coUisional 
properties closely constrain important features of the disk. For the more dilute 
disks, H is directly correlated with c and seems to adhere quite well to the 
dilute estimate oi H k, c/VL. In contrast, it behaves quite differently in the 
case of Vc/v}, — 1. Here H possesses a subtle turning point at a critical r 
after which it begins to increase. The increase is caused by the vertical fiux 
of momentum due to collisions, the 'collisional pressure', which increases with 
r. As the translational pressure drops and the disk vertically contracts this 
contribution will become important and will work to fatten the disk. If we 
had included self-gravity, however, the collisional pressure would be met and 
overcome and the disk would flatten further on increasing r (see AT86). At 
some point, though, this compression must stop. 

Fig. 3c is interesting because it shows the disk approach the dilute energy 
balance as Vc is increased. The limiting £(t) curve in this case being that of 
Goldreich and Tremaine (1978) (cf. their Fig. 2). In the dense regime (low Vc) 

the injection of energy by the collisional stress source becomes important and 
so each collision must dissipate more energy than otherwise. This explains the 
low values of e the energy balance requires when Vc/vi, = 1 and 0.25. . 

Figs 4a and 4b present the nonlocal and local kinematic viscosity against r. 
Here we find relatively good agreement with Salo (1991) for each Vc/vb- As with 
the velocity dispersion, the kinetic treatment overestimates the viscosities, 
especially the nonlocal viscosity for larger ratios of Vc/vi,. Otherwise, the figure 
displays quite well, in the relative 'convexity' and 'concavity' of the curves. 
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Fig. 3. Various equilibrium properties computed by a polytropic kinetic model 
(a = 0.9) with a preaveraged e (solid curves), and by the simulations of Salo (1991) 
which we denote by 'S91' (triangles). We employ different ratios of Vc/v^ (inserted 
adjacent to the relevant curve) but fixed p = 0.234. The quantities are plotted as 
functions of optical depth r. The four panels correspond to: (a) equilibrium velocity 
dispersion C, (b) disk scale height H, (c) averaged equilibrium coefficient of restitu- 
tion e, and (d) central filling factor FFq. The Goldreich and Tremaine (1978) dilute 
disk energy balance is referred to as GT78 in subplot (c). 

the different behaviours associated with the two regimes of disk — the warm 
and dilute on one hand, and the dense on the other. 

In Fig 5 we plot the quantity ru, which is proportional to the angular mo- 
mentum flux and hence a key indicator of the stability of the equilibrium. We 
present three interesting cases. The first, corresponding to Vc/vb = 20, repre- 
sents a 'bimodal' profile; such cases have been discussed by Hameen Antilla 
(1982) and Spahn and Schmidt (2006), and were first revealed in simulations 
by Salo (1991). For low to middle r the disk is in the dilute regime; there 
the angular momentum flux increases sharply and then gently decreases after 
reaching a turning point. As explored in Latter and Ogilvie (2006), the region 
of large positive slope does not lead to viscous overstability because the stress 
(dominated by the translational component) is nonlocal in time; however, the 
region of negative slope should give rise to the viscous instability. At high 
optical depths the disk enters the dense regime and the collision frequency be- 
comes sufficiently large to engender a significant collisional stress: the angular 
momentum flux begins to increase gradually under its influence. 
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Kinetics 

■< S91 




Fig. 4. Local and nonlocal components of the kinematic viscosity. Panel (a) presents 
the nonlocal component of the kinematic viscosity v^^, and (b) the local compo- 
nent of kinematic viscosity v^. As earlier different ratios Vc/vh are employed with 
p = 0.234, and the solid curves represent the polytropic kinetics and the triangles 
Salo's simulations. Both viscosities are scaled by a^Q. 

40 I 1 , , , 1 
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Fig. 5. The total angular momentum flux tu = t{u^^ + v^) versus r at different 
ratios of Vc/vh- The flux is scaled by Q.a?. 

The second curve corresponds to fc/fb = 10, and it shows the coUisional stress 
asserting itself at a lower r, effectively negating the 'window' of decreasing 
angular momentum flux exhibited by the Vc/vb = 20 curve. In this intermediate 
range the slope of the curve is near to zero as the opposite tendencies of the 
local and nonlocal stresses cancel. The disk is perhaps marginally stable to 
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Fig. 6. The scaled collision frequency Wc/f^ as a function of optical depth r for two 
illustrative ratios of Wc/'ffe- 
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P 

Fig. 7. The equilibrium velocity dispersion at FFq = 0.2 as a function of the expo- 
nent p with fixed Vc = V},. 

the viscous instability in this range. 

The third curve corresponds to a disk which lies in the dense regime for all 
T. Here u is always dominated by the nonlocal component. In this case the 
angular momentum flux gradually increases and may at a critical r possess a 
gradient sufficiently large to initiate the viscous overstability. 

Fig. 6 plots the collision frequencies for two ratios, Vc/vb = 40 and 1. These 
illustrate the dilute regime and the dense regime respectively. In the latter 
case Uc clearly exhibits its linear dependence on r, while the former shows 
a gentle superlinear dependence for sufficiently large r, this mainly arising 
from the Enskog factor. The kinetic model uses expression (62) to evaluate 
the collision frequency. The matching between the two approaches is good in 
general, though at larger r, there is an increasing discrepancy in the dense 
case. 

Finally, in Fig. 7 we plot the variation of c with respect to the elasticity 
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exponent p for fixed Vc — Vf, and FFq = 0.2. The trend reveals that increasing p 
decreases the equihbrium velocity dispersion. This follows because an elasticity 
law with a larger p varies with impact velocity more steeply. With such a law, 
a collision characterised by a given impact velocity is more dissipative than if 
the slope were more shallow. If collisions are more dissipative in general then 
the equilibrium c may be smaller both to minimise the collision frequency and 
to maximise the injection of energy by the collisional stress. 

3. 9 Importance of vertical nonlocality 

Having applied the kinetic model to equilibrium calculations, compared the 
properties computed with those from simulations, and discussed these prop- 
erties more generally, we shall now test some of the assumptions it uses. First, 
we examine the approximation of vertical locality, introduced in Section 2.2.6. 

In order to best test the assumption we present the equilibrium characteristics 
of two representative disks: one without self-gravity and one with self-gravity. 
These two cases correspond to relatively thick and thin disks, respectively, and 
hence the assumption of vertical locality will impact on each quite differently. 
Otherwise, the model incorporates 

• a preaveragcd e, the piecewise power law, using (45) and (46) and with 
p = 0.234 and Rc = 2 

• Gaussian vertical structure. 

We employ Gaussian vertical structure because this model admits analytical 
vertical averages in the vertically nonlocal case, in contrast to the polytropic 
model (see Section 3.4). We believe the results we obtain with the Gaussian 
profile would also hold qualitatively in the polytropic case. Recall that the 
only change made by the assumption of vertical locality is omission of the 
'nonlocality factor', exp[—(a'^/H'^)kl] in the collision integrals. 

First we set Gg = and compute non-self-gravitating equilibria with and with- 
out vertical nonlocality. Selected equihbrium properties, as functions of FFq, 
are plotted in Fig. 8. The two models agree closely because of the relatively 
large scale height, H, plotted in Fig. 8b. This quantity is approximately 2.5 
at its smallest; consequently, the 'nonlocality factor' varies at most between 
about 1 and 0.85, and much less at higher central filling factors. The discrep- 
ancy between the two models can then be at most some 15%, and it appears 
its effect is much less. 

Next we set Gs = 6.6, a realistic value for solid ice paticles in the B-ring. The 
equilibrium properties are plotted in Fig. 9, as before. The discrepancies be- 
tween the two methods are more pronounced, and especially attending those 
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quantities most closely aligned with the vertical structure, r and H. These two 
deviate substantially at large central filling factors as then the disk semithick- 
ness approaches a particle radius. For these very thin disks the 'nonlocality 
factor' can vary between 0.367 and 1 and thus play a significant role. However, 
we are encouraged by the agreement exhibited by the horizontal properties, 
such as velocity dispersion, which are not nearly as affected. 

The reason why the vertical properties are more affected than the horizontal 
ones is because of the local model's overestimation of the vertical coUisional 
pressure P^^ in the height equation, which is were H and hence r are deter- 
mined effectively. In a very thin disk P^^ will be very small: most colliding 
particles are in the same plane and hence nearly all collisional momentum 
transport will lie in this plane. But the collisional dynamics of the local model 
is insensitive to the vertical stratification and hence to this fact. Evidence for 
this interpretation can be observed in Fig. 9b where we see the disk semithick- 
ness H increase at a lower FFq than it should because of the overestimation 
of the vertical momentum fiux. 




Fig. 8. Comparison between the vertically local model (the solid line) against the 
vertically nonlocal model (triangles) in a non-self-gravitating disk, Gs = 0. In (a) 
we plot selected 'horizontal properties' of the equilbria (C, P^y, Hxy), and in (b) the 
'vertical properties' {H, r). These are presented as functions of central filling factor 
FFq. Otherwise we have adopted Gaussian vertical structure and a preaveraged e 
with p = 0.234 and Rc = 2. 

In summary, however, the assumption of vertical locality is certainly justified. 
Very little error is introduced in disks which possess a H > 2a, and the hori- 
zontal properties are approximated adequately in self-gravitating disks where 
the semithickness approaches a particle radius. It is only the optical depth r 
and H itself in very thin disks that we must be careful interpreting when FFq 
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Fig. 9. As for Fig. 8 but for a self-gravitating disk, Gg = 6.6. 

is high. This is pleasing because the assumptiom of vertical locality provides 
a number of mathematical advantages. Particularly useful is its decoupling 
of the height equation from the horizontal pressure tensor equations and the 
simplification of the integrands. These simplifications, also, are crucial in the 
polytropic model, and also in the interpolation, and expansion models, as they 
facilitate the complete integration of the collision terms. 

3.10 Validity of preaveraged coefficient of restitution 

In this section we test the preaveraging approximation against a model which 
incorporates the variation of e consistently within the collision integrals. The 
coefficient of restitution is set as a straight power law: e = (g ■ k/vc)~^ where 
Vc = 0.01 cm s~^ and p = 0.234. This functional form supplies analytical 
expressions for both {e)q and Aj{e). Expressions for the former are listed in 
(42) and (43) and the latter are simple power laws but with a slightly different 
'Wc' (we omit these expressions). Thus our choice of the straight power law 
ensures that the integrations we need to undertake are only two-dimensional 
and not three, as would be the case for the piecewise function. In addition we 
set the velocity scale aQ = 0.02 cm s~^ and Gg = 0. Otherwise, the model 
assumptions are 

• polytropic vertical structure with a = 0.9 

• vertical locality 

In Fig. 10 we plot selected equilibrium properties. The agreement is quite 
good generally and excellent for the coUisional stress. The preaveraged model 
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Fig. 10. Selected equilibrium properties as computed by a model which incorporates 
a preaveraged e and a model that allows e to vary within the collision integrals. A 
straight power law was adopted for e, with p = 0.234 and Rc = 2. The quantities 
H and C appear in panel (a), and the collisional and translational stresses P^y and 
lixy in panel (b). The vertical structure is assumed to be polytropic. 

slightly overestimates the equilibrium velocity dispersion and collisional pres- 
sure tensor. Overall, however, these results reveal that this approximation 
gives excellent qualitative and good quantitative agreement. 

3.11 Interpolation model 

Finally, we investigate the interpolation model introduced in Section 2.2.8 and 
worked out in detail in Appendix A. The F functions were interpolated on 
the range ^ G [—0.75, 0.75] with a mini-max method which returned a relative 
error of at most 3%. The Fij which appear in Eqs (47) were computed as 

Fii = 1.9619, Fi2 = 5.3174, F13 = 6.545, Fu = 3.5450, 
F21 = 1.0000, F22 = 1.7725, F23 = 0.9352, 
F31 = 1.7725, F32 = 4.5401, F33 = 3.5450 

(to 4 decimal places). Also we assume 

• polytropic vertical structure with a = 0.9, 
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Fig. 11. Selected equilibrium properties as functions of r computed by the full model 
(points) and the interpolation model (solid line). Both models adopt polytropic 
vertical structure (a = 0.9) and a preaveraged e with p = 0.234 and Rc = 2. 
The velocity dispersion C is plotted in panel (a), the semithickness H in (b), the 
coUisional stress P^y in (c), and components of the translational stress 11 in (d) 

• vertical locality, 

• a piecewise power law, which is preaveraged with p = 0.234 and Rc = 2 

Because the collision integrals are linear in a the equilibrium problem is greatly 
simplified. We can solve Eqs (48)-(51) analytically for the xx, xy and yy 
components of a, leaving a single nonlinear algebraic equation for c in terms 
of the governing parameters, r, Rc, p, Gs- The equilibrium calculation is then 
rapidly completed by a one-dimensional Newton-Raphson algorithm. 

A comparison of the equilibrium properties as computed by the interpola- 
tion and exact method are presented in Fig. 11. The agreement here is good 
throughout the full range of r which is somewhat surprising: we expected dis- 
crepancies near r = arising from the truncation in anisotropy (a). Perhaps 
this accounts for the deviation of disk semithickness H at low r, but it is puz- 
zling that this vertical property is singled out. In general, given the truncation 
in a, the excellent agreement for general r suggests that dense rings exhibit a 
level of anisotropy which can be captured adequately by a linear model. 

Because this method reduces the collision integrals to algebraic expressions 
(albeit complicated), more advanced linear and nonlinear analyses of the ring 
kinetics are easier to undertake. Moreover, its success in equilibrium calcula- 
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tions encourages its application in these investigations. 



4 Linear stability 

Now that the statistics of the dense disk equihbrium have been estabhshed, we 
shall examine its response to small perturbations. Our interest is principally 

in the growth rates and stability of the overstablc modes as functions of the 
kinetic parameters. We will not supply an exhaustive survey in this paper but 
rather illustrate the principal features while also touching on a comparison 
with the simulations of Salo et al. (2001). 



Linearised equations 

The disk admits steady homogeneous solutions characterised by the constant 
optical depth tq and velocity dispersion tensor W°. Suppose that we perturb 
these equilibria with an infinitesimal axisymmetric disturbance 



The system is subsequently nondimensionaliscd along the following lines: time 
with respect to Q^^, space with respect to Co/f2, velocity with Co, and velocity 
dispersion tensor with Cq. We now use the nondimensionaliscd quantities but 
with no change in notation. 

The linearised equations read 



W 



U 



r 



To + T'{x,t), 

W° + W'(a;,t), 
$°, + $'^(a;,t), 
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where the perturbed coUision terms, Q'^^ and (Pij)', are expressed in the form 

\ Pirr ] ^ \ I '''' \dW^y ) ^ 




dW, 





and Cij is the rate of strain tensor. To these we add the perturbed Poisson 
equation, in order to obtain and the perturbed height equation, in order 
to relate the perturbed central filling factor FFq with r'. 

Once we express the radial perturbation as Fourier modes oc exp{ikx + st), 
we let the nondimensionalised solution of the linearised Poisson equation take 
the form 

^'d - -^r' (70) 
(Binney and Tremaine, 1987), where 

and is equal to (Qtq)"^, where Q is the Toomre parameter. Once this is substi- 
tuted we obtain a seventh order algebraic eigenvalue problem for the growth 
rate s. This in turn supplies a seventh order dispersion relation of the same 
structure as Eq. (25) in Latter and Ogilvie (2006), but in which the coeffi- 
cients are functions of Gg, g, k, the equilibrium values of W and P"^, and the 
derivatives of W and P'^. If we adopt the general dense formahsm devised in 
Section 2.2.6, evaluating these coefficients requires the numerical calculation 
of a number of surface integrals and their derivatives. 

Before continuing, we will say a few words about the parameter g given by 
Eq. (71). Though it appears that g simply derives from the quantities Gg and 
R, a quick calculation shows that it can take large values in a typical dense 
equilibrium, and hence Q can take very low values. For example, if Gs = 6.6, 
a = 1 m, r ~ 1 and G/{aQ) ~ 1.5 (see Fig. 9), we compute Q ~ 0.23. Kinetic 
models characterised by these low values of Q exhibit viscous overstability 
on arbitrarily short scales. Some aspects of this issue have been examined 
in Latter and Ogilvie (2006) in the context of viscous instability. Part of 
the problem originates in the small velocity dispersion of axisymmetric self- 
gravitating equilibria, an admittedly artificial scenario. A real self-gravitating 
disk will be subject to non-axisymmetric instability for a Q less than about 2, 
and the recurrent wake structures it engenders will heat the disk sufficiently to 
keep Q above unity. In addition, gravitational scattering and focusing (effects 
we neglect) will supply extra heat sources. Self-gravitating simulations bear 
out these point and show that Q is typically unity or higher (Salo, 1992a). 
The second important issue is our omission of the third order moments which 



43 



are important on short scales and which should extinguish instability. These 
considerations will not be addressed in this paper, but we will revisit it in a 
future nonlinear axisymmetric study where a solution to this problem is more 
urgent. 



4.2 Results 



The model we employ incorporates 



• poly tropic vertical structure with a = 0.9, 

• vertical locality, and 

• a preaveraged piecewise power law for e. 

Both the equilibrium and its linear stability depend on six parameters: Tq, p, a, 
m, Q, and Vc, all of which can be determined by observation or experiment. The 
nondimensionalised system is parametrised by the dimensionless quantities r, 
p, Rc and Gg- 

The perturbed collision terms depends on the perturbations of FFq via the 
Enskog factor. But this must be related to the other perturbations before we 
can proceed. The relevant equation can be obtained via (58): 
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This equation, which relates the perturbations of FFq, H and r, conveys in- 
formation pertaining to vertical disk structure into the linear dynamics via 
the equilibrium H and its derivatives. In particular, it introduces vertical self- 
gravity into the problem. Self-gravity also appears in the x component of the 
momentum equation, (64). 

Once the linearised equations are sorted, our procedure is to step through 
FFq, evaluate the equilibrium at each point and obtain the dispersion relation, 
which is then solved for various k (or wavelength). We perform three runs: 

• without self-gravity in either the height equation or in the momentum equa- 
tion, Fig. 12; 

• with self-gravity in the height equation alone (thereby modelling only its 
enhancement of the vertical restoring forces). Fig. 13; 

• with self-gravity in both equations. Fig. 14. 
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The algorithm evaluates the growth rates of the seven linear modes discussed 
in Latter and Ogilvie (2006). These are: the two potential viscous overstable 
modes, the potential viscous instability mode, the energy mode, the two po- 
tential anisotropic overstable modes, and the relaxation mode. We, however, 
concentrate on the two potentially overstable ones (left and right-going waves) . 
First, their curves of marginal stability are plotted in the (r. A) plane, where 
A = 27r/A; is wavelength. Then, to these we add the curves of constant growth 
or decay, which we truncate as they collapse onto the line of marginal stability. 

We adopt the Bridges et al. (1984) elasticity law, so as to compare with the 
simulations of Salo et al. (2001); thus p = 0.234 and Vc = 0.0077 cm s~^. 
The particles are assumed to possess radii of 100 cm and the shearing sheet is 
located at a radius of 100, 000 km. As a result afl = 0.0195 cm s^^ and Rc = 
afl/vc — 2.53. We vary Gg so as to compare our results with the simulations, 
in which the particles' internal density took several values less than 900 kg m~^ 
(solid ice). In addition, Salo et al. conducted a number of simulations where 
the vertical component of self-gravity was mimicked by an enhancement of the 
particles' vertical epicyclic frequency, flz- The effective Gs associated with such 
as enhancement canot be readily determined, though we have the estimate 

+ {dl<PD)z=o + 3Gs FFo), (72) 

in which it has been assumed that (po varies much slower with x than with z. 
Typically, Salo et a/.'s simulations set Qz/^ = 3.6 which corresponds to very 
large Gg, especially for low FFq. When FFq is about 0.3, Gg is still greater 
than 10. This is surely an overestimate, and Eq. (72) should be treated with 
some caution. 

4-2.1 Growth rates 

We find that a non-self-gravitating disk is stable to the viscous instabihty for 
all the optical depths we examine: < r < 5. This conforms with the dense 
gas A^-body simulations. Above a critical r the viscous overstability occurs. 
This value is approximately r = 2.5 which corresponds to a central filling 
factor of 0.375 (cf. Fig. 12). In contrast, the A^-body simulations of Salo find 
overstability only for very high optical depths, of order 4. This disagreement is 
not unexpected: high filling factor disks in monodiperse simulations can 'relax' 
their close packing via the formation of layers or 'sheets'. Such disks, as a 
consequence, can possess filling factors below the critical level of overstability, 
despite being quite optically thick. 

The stability curves of the long modes differ very little from each other: once 
the critical r is reached nearly all the overstable modes become unstable, 
which mirrors the hydrodynamical stability curves of Schmidt et al. (2001). 
As the disk makes the transition to instability, modes on lengths of some 200 
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m will grow fastest. The shortest overstable scales are approximately 80 — 90 
m. These are a little longer than usually predicted by hydrodynamics and 
simulations. 
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Fig. 12. The curves of constant growth and decay rates, Re(s), for the overstable 
modes in the (r, A) plane in a non-self-gravitating disk (Gs = 0, g = 0). The 
model adopts polytropic vertical structure and a pre- averaged e with p = 0.234 and 
Rc = 2.53. Negative growth rates are plotted with dashed lines and positive growth 
rates with solid lines. The curve of marginal stability, Re(s) = 0, is the dotted line. 



In order to examine the contribution of the vertical component of self-gravity 
we set Gs = 3.0 in the height equation, but zero in the momentum equation. 
This captures the vertical compression caused by self-gravity (and the subse- 
quent increase in collision frequency and its rate of change) but will not alter 
the structure of the dispersion relation. The value adopted for Gs corresponds 
to internal particle densities of pp = 408kgm~^, which is less than solid ice 
but is illustrative nonetheless. 

The curves of constant growth and decay are plotted in Fig. 13. These are 
similar to those in Fig. 12; the chief difference lies in a shift to lower optical 
depths and slightly shorter scales. The critical r in this case is ~ 1.15 and the 
critical FFq approximately 0.335. 

The fastest growing modes at the stability transition are now on scales of about 
150 — 200 m, and the shortest affected scales are about 60 — 70 m. These agree 
a little closer with the generic hydrodynamic profiles. 
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Fig. 13. The curves of constant growth and decay rates Re(s) for the overstable 
modes in the (r, A) plane, in a disk where self-gravity is included in the height 
equation but not elsewhere. We have set = 3 and g = 0. 

Next we include self-gravity in both the height equation and in the x compo- 
nent of the momentum equation. We seek to compare directly with the self- 
gravitating run of Salo et al. in which pp = 225 kg m~'^. Therefore, Gg = 1-65 
by Eq. (56) and g takes the appropriate value from Eq. (71). The subsequent 
curves of constant growth and decay are plotted in Fig. 14. 

With the addition of self-gravity in the momentum equation it is the interme- 
diate lengthscales (~ 100 m) which become unstable first, as illustrated by the 
pronounced salient in the curve of marginal stability. The longest wavelengths 
are unaffected. On the other hand, decreasing Gs from 3 to 1.65 increases the 
critical r of overstability on the longest scales, from about 1.15 to approxi- 
mately 1.43 (which corresponds to a central filling factor of 0.34). In contrast, 
the lowest critical r for overstability on the intermediate scales is about 1.23 
(for A ^ 100 m). 

These results are in good agreement with the simulations, which find that 
overstability begins at approximately 1.2 (see Fig 5a in Salo et al, 2001) and 
on scales of 100 — 150 m. This conformity contrasts with the poor agreement 
in the non-self-gravitating run and proceeds essentially from the introduction 
of the vertical self-gravity. In simulations the self-gravity negates, to some 
extent, the relaxation of packing fraction and allows the higher filling factors 
necessary for overstability. It thus removes the obstacle that prevented the 
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Fig. 14. The curves of constant growth rate Re(s) for the overstable modes in the 
(r, A) plane, in the case when self-gravity is properly included in both height and 
momentum equations. Here Gg = 1.65 and g = GsR, cf. Eq. (71). 

kinetic model and the simulations agreeing. 

Finally we set Gg = 6.6, corresponding to solid ice particles, Pp = 900kgm~^. 
The parameter g was set to its appropriate value from (71), and thus the 
Toomre parameter was generally very small. The critical r is now shifted 
to about 0.8 for long scales and somewhat less for intermediate scales. As 
mentioned earlier, in Section 4.1, we find that the overstability is not quenched 
on the short scales, thus the 'bulge' in the stability curve associated with 
self-gravity (cf. Fig. 14) extends down to lengths of a particle radius. This 
behaviour is surely unphysical and derives not only from the small Q but also 
the omission of the third order moments in the pressure tensor equation. 

4-2.2 Marginal curves 

Finally we let the elastic properties of the particles vary by letting p and the 
parameter Rc = aQ/vc take other values. Marginal curves of Re(s) = for 
viscous overstability are computed in {p, r) space for given k and Rc for a 
self-gravitating disk with Gg = 6.6 but with g = 0. These curves are obtained 
from a simple (but robust) bisection root-finding method. Fig. 15 plots three 
illustrative curves for a mode with k = 0.01, which corresponds roughly to 
A = 600 m. The parameter Rc takes the values 2.53, 0.67, and 0.25. The first 
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two correspond to 100 cm size particles at a radius of 100, 000 km in the B- 
ring with the elasticity data of Bridges et al. (1984) and Supulver et al. (1995) 
respectively. 




Fig. 15. Marginal curves of overstability in the (p, r) plane for k = 0.01, Gs = 6.6 
and various Rc = aVi/vc- Regions above the curves are overstable, regions below are 
stable. 

The stability trend here is easy to interpret: the larger r, p, and Rc the more 
vulnerable the disk to viscous overstability. The last correlation derives from 
the fact that the greater Rc the more important the background shear flow, and 
consequently the greater the collisional viscous stress. As Rc — >■ the nonlocal 
contributions vanish. The effect of p lies in its control of the equilibrium disk 
velocity dispersion. The larger p the lower C, and hence the more important 
the collisional viscous stress (cf. Fig. 7). 

Another fact to note is the sensitivity of the critical r to the elasticity laws. 
The Bridges et al. data (for which p = 0.234) gives a critical r of roughly 0.8, 
but the Supulver et al. data (for which p = 0.19) gives a value nearer to 1.5. 
In contrast, in these two cases and more generally, the critical FFq does not 
change as dramatically and lies somewhere in the interval 0.32 — 0.375. We 
conclude that the viscous overstability responds essentially to a sufficiently 
large filling factor, the critical value of which does not vary significantly with 
respect to the collision law (or self-gravity) . It rather than r is the key param- 
eter governing the onset of instability. 

Though the stability results presented in Fig. 15 hold only for a mode of length 
some 600 m, the stability of all the scales should follow qualitatively the same 
trend, irrespective of the precise value of Gg and g. Indeed these results are 
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consistent with the self-gravitating simulations of Salo et al. (2001) when they 
vary (cf. their Fig. 5b). 

Given that optical depth is a key observational property, and that it depends 
so delicately on the elastic properties of the particles, our detailed model is 
a useful first step in constraining the collisional parameters of ring particles. 
Certainly the model needs refining, particularly with respect to self-gravity, 
but it does provide an account of dynamical processes which can be compared 
directly with observations. Studies such as Salo and Karjalainen's (2003) com- 
bined dynamical and photometric modelling have already made inroads into 
such comparisons. 



5 Conclusion 

In summary, we have developed a workable formalism with which to attack the 
dynamics of a dense ring, one in which the collisional transport/production 
processes and filling factor effects are fully included, as is the dependence of e 
on impact velocity. The model both simplifies and generalises the work of AT86 
and as a consequence encourages us to conduct more advanced dynamical 
analyses. 

The theory is founded on that of Enskog kinetics and thus shares its assump- 
tions regarding the neglect of ternary or higher collisions and the modelling 
of space and velocity correlations. An additional assumption is the limiting 
of the velocity dispersion's spatial variation to scales significantly longer than 
a particle radius (i.e. vertical isothermality) . Also the velocity distribution is 
assumed to be a triaxial Gaussian, and the vertical structure of the number 
density to take a predetermined form, either a Gaussian or a polytrope. The 
theory must be supplemented, with the functional forms of the Enskog factor 
and the coefficient of restitution. 

The full model can be further simplified if additional assumptions are made, 
such as those of locality in density, a preaveraged coefficient of restitution, or if 
the collision integrands are interpolated and anisotropy assumed small. Each 
of these models fared well in comparisons with the exact approach (Sections 
3.9 — 3.11). Thus we employed them generally. 

We undertook a comparison of equilibria computed by the kinetics of Araki 
and Tremaine (1986) and the simulations of Wisdom and Tremaine (1988) 

and Salo (1991). The former two comparisons exposed the shortcomings of 
the Gaussian model of vertical structure at moderate and high central fill- 
ing factors in monodisperse disks, and revealed an error in the numerical 
implementation of the Araki and Tremaine model. It was also shown that a 
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polytropic model of vertical structure performed well in these regimes and 
acceptably at lower filling factors. 

The polytropic kinetic model was tested in equilibrium calculations against 
the results of Salo's (1991) A^-body simulations. Overall the agreement was 
good, with the major discrepancies within the range of error introduced by 
the various aproximations. 

Next we probed the linear theory of the kinetic model and reproduced quali- 
tatively much of the behaviour presented by the particle simulations of Salo et 
al. (2001). We find that without self-gravity a disk composed of 100 cm radius 
particles exhibits viscous overstability at a critical optical depth of about 2.5. 
This is a significantly lower value than Salo et aVs, which is approximately 4, 
and proceeds essentially from the kinetic model's insistence on a fixed vertical 
structure. A simulation can 'relax' its close packing by the formation of vertical 
sheets and thus remain below the critical filling factor necessary to instigate 
overstability. (Note however that the sheets are artefacts of a monodisperse 
model and hence unrealistic.) This discrepancy vanishes when self- gravity is 
properly included and agreement between theory and simulations is quite good 
for comparatively light particles. In this case modes of 100 m are the most un- 
stable and possess a critical optical depth of about 1.23. In general, inclusion 
of the vertical component of self-gravity decreases the critical optical depth 
above which overstability manifests, while inclusion of the radial component 
of self-gravity induces overstability preferentially on intermediate scales, while 
leaving the longest scales untouched. 

When the self-gravity is large, for instance when particles are as dense as solid 
ice, the linear dynamics on short scales is incorrectly described. The viscous 
overstability then extends to lengths of about a particle radius, which is un- 
physical. This arises because of the omission of third order moments from the 
pressure tensor equations and also because of the very low velocity disper- 
sions which characterise axisymmetric disks in which gravitational collisions 
have been neglected. In reality, non-axisymmetric wakes, and gravitational 
collisions, warm the disk sufficiently and on short scales 'thermal diffusion' 
extinguishes instability. This issue must be addressed more fully in future 
nonlinear studies which employ the kinetic model and assume axisymmetry. 

Within the model assumptions, we predict that overstability can occur at any 
radius of Saturn's rings provided the optical depth exceeds a critical value. 
This value depends rather sensitively on the parameters of the collision law, in 
our case p and v^, in the B-ring, a plausible range of critical r is between 0.8 and 
1.6 (cf. Fig. 15). It must be said, however, that these estimates omit the effects 
of the gravitational wakes and size distribution both of which militate against 
the overstability's onset. Even so, the dense B-ring seems the most likely 
place in which overstability occurs, as it exhibits unambiguously high optical 
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thicknesses. Finescalc structure in such optically thick regions we conclude 
should be associated with this instability. 



Our kinetic model could be developed in a number of ways. For instance, it 
would be straightforward, though tedious, to include the spin degrees of free- 
dom and, possibly, size distribution. In the former case, we may use similar 
techniques to simplify the spin collision integrals. If a similarly workable for- 
malism was developed we could examine, amongst a great many things, the 
action of particle spin on the onset of the overstability. Simulations have al- 
ready begun probing this problem and await a kinetic comparison (Morishima 
and Salo, 2006). However, it should be noted that at present the frictional 
properties of ring particles are poorly understood. Also, the effect of gravi- 
tational encounters between particles should also be included. This is an im- 
portant but difficult task. At present the model omits this important heating 
mechanism and thus underestimates the velocity dispersion of self-gravitating 
equilibria. 



There are a numerous ring-related problems on which the full or interpolation 
formalism could be put to use immediately. The most obvious is the nonlin- 
ear saturation of the viscous overstability. The seven evolutionary equations 
could be numerically simulated to investigate the long-term evolution of the 
instability and how the elasticity parameters influence the statistics of the 
quasi-steady state to which it leads. The damping of nonlinear density waves 
launched by moons such as Prometheus, Pandora, etc could also be studied, as 
well as the wavy edges and wakes induced by Daphnis and Pan in the Keeler 
and Encke gaps. Theoretical models and their predictions have suffered from 
an approximate description of the viscosity in the ring. The full dense gas 
model could supply a more correct account of these effects. These would be 
formidable computational problems as at each time step and each spatial grid 
point we would need to evaluate a number of two-dimensional integrals. Such 
projects may best be tackled with the interpolation model with its algebraic 
expressions for the collision terms. That said continuum models, such as those 
we develop here, can reach the longer length and time scales necessary to fully 
describe the nonlinear behaviour. A^-body methods are still too computation- 
ally intensive, particularly when self-gravity is added, to feasibly capture these 
scales. 
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6 Appendix: Interpolated collision terms 



Presented here are the collision terms computed from the interpolation method 
to order a. We have: 



L2 = Anc^ jfFao (I2 + fa) - — F31 (ae + ea + 2e + IV(e)a + IV(e)l2) 
2 

-— ^F32 (8[e^aL + ae'] + 4Tr(e) [ea + ae] + 2[Tr(e)2 + 2Tr(e')]a 
+30l6:e'-fl8:a:e')}. 



The Fij are pure numbers which are determined from the interpolation method 
and interpolation range; I2 is the rank 2 identity tensor and the other tensors 



q = - a'(l + e)Y [(1 + e)Li - 2L2] 



(73) 



where the Li tensor is defined by 




and the L2 tensor by 
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are defined through 
I4 : a = |a, 

I4 : e = |e + iTr(e) I2, 

le : a : e = ^ {2(ae + ea) + IV(e)a + 'IV(ea)l2} , 

le : e' - {Se^ + 4Tr(e)e + Tr(e)'l2 + 2Tr(e2)l2} , 
Ig : a : e^ = |8(ae^ + eae + e^a) + 4Tr(e)(ae + ea) 

+4Tr(ae)e + [Tr(e)2 + 2Tr(e2)]a + 2Tr(e)Tr(ae)l2 + 4Tr(eae)l2} 
Is : = ^ {48e^ + 24Tr(e)e2 + 6[Tr(e)2 + 2Tr(e^)]e 

+[Tr(e)3 + 8Tr(e^) + 6Tr(e)Tr(e')]l2} . 

The coUision momentum flux tensor is 

= 80F(1 + e)Y a' {F20 {^h + : a) 

-^^-(|l4:e+^Ia:a:e) + ^F22la:e^}. 

If the Fjn are set to zero except for Fio = 2, F20 = ^/7l, and F30 = 1 then we 
recover the dilute gas formahsm of Goldreich and Tremaine (1978) to order a. 
See Section 2.2.9. 
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